Fundamental high pressure calibration from all-electron quantum Monte Carlo calculations 
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We develop an all-electron quantum Monte Carlo (QMC) method for solids that does not rely on pseudopo- 
tentials, and use it to construct a primary ultra-high pressure calibration based the equation of state of cubic 
boron nitride(c-BN). We compute the static contribution to the free energy with QMC, and obtain the phonon 
contribution from density functional theory, yielding a high-accuracy calibration up to 900 GPa usable directly 
in experiment. Furthermore, we compute the anharmonic Raman frequency shift with QMC as a function of 
pressure and temperature, allowing optical pressure calibration in table-top experiments. In contrast to present 
experimental approaches, small systematic errors in the theoretical EOS do not increase with pressure, and no 
extrapolation is needed. This all-electron methodology is generally applicable to first-row solids, and can be 
used to provide a new reference for ab initio calculations of solids and to benchmark pseudopotential accuracy. 
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Although the number of studies of materials using density 
functional theory (DFT) continues to grow explosively, the 
accuracy of their predictions is variable, sometimes excellent 
and sometimes poor, limiting the confidence one may place 
in DFT results as a quantitative calibration for experimen- 
tal studies. Quantum Monte Carlo (QMC), particularly dif- 
fusion Monte Carlo (DMC), is the highest- accuracy method 
for finding the ground state of a many-electron Hamiltonian. 
For solids with atoms heavier than He, however, the Hamilto- 
nian itself is approximated using pseudopotentials (PPs) based 
on a lower-accuracy theory, limiting its reliability, and as we 
show, commonly used PPs give disparate results. In this Let- 
ter, we push the state of the art in accuracy by introducing a 
method for an all-electron DMC simulation of solids, elimi- 
nating the bias from pseudopotentials, and apply it to create a 
high-accuracy pressure calibration scale. 

The combination of ultra-high pressure mineralogy with 
seismology has yielded a wealth of insight into the internal 
structure of our planet. Pressure is the key that links these 
disciplines, mapping phase transitions and mineral properties 
to planetary depth. Establishing an absolute pressure cali- 
bration at multi-megabar pressures, however, poses a funda- 
mental and continuing problem for high-pressure experiment. 
Current primary calibrations are based on data from shock- 
wave experiments, which infer pressure from conservation of 
momentum and energy as the shock traverses the sample. Ex- 
perimental uncertainty, combined with errors in model extrap- 
olation, yield scales which differ from each other by as much 
to 7% at room temperature, with even greater discrepancies at 
high temperature[l]. Such disparity remains a serious obsta- 
cle to a quantitative understanding of Earth's interior. 

A pressure calibrant is a material with a known equation of 
state (EOS), a small sample of which may be included in hy- 
drostatic equilibrium with the test subject (e.g. ruby[2]). As 
an alternative to shock experiments, high-pressure Brillouin 
scattering, which provides a measurement of the bulk mod- 



ulus, in conjunction with X-ray diffraction measurements of 
volume, can be integrated to provide an EOS, but a correc- 
tion must be made to transfer from an adiabatic to an isother- 
mal path[3]. New approaches such as quasi-adiabatic Z-pinch 
based experiments [1, S\ also hold future promise for a pri- 
mary scale. There have been attempts to refine the ruby 
scale|0|], and new calibrations have also been suggested ||3fl. 
Cubic boron nitride(c-BN) has been identified as a promising 
material for a new scale ||7|]. In this Letter, we provide a new 
pressure scale based on DMC calculations of the EOS and Ra- 
man frequency of c-BN. This theoretical approach has the ad- 
vantage that, in contrast to present experimental approaches, 
the method works equally well under high compression, and 
uncertainty does not grow with pressure. 

Pressure is the negative volume derivative of the Helmholtz 
free energy. In a wide-gap insulator such as c-BN, the free 
energy can be written as a sum of the frozen lattice enthalpy, 
dependent only on volume, and a phonon thermal free energy, 
which depends on both volume and temperature. Since the 
static enthalpy is the dominant contribution at ordinary tem- 
peratures, errors in a theoretical EOS can most often be at- 
tributed to the static part. Previous calculations of the EOS 
of c-BN have been based on DFT[7], which use approximate 
functionals to treat electron exchange and correlation. Several 
functionals are in common use, each giving rise to a different 
EOS, and there is no a priori way to predict which will give 
the most reliable result. 

Quantum Monte Carlo simulation methods explicitly treat 
electron exchange and correlation instead of resorting to ap- 
proximate functionals. Variational Monte Carlo (VMC) com- 
putes properties by a Metropolis sampling of a trial wave func- 
tion. Diffusion Monte Carlo samples the many-body ground 
state of the Hamiltonian through a stochastic projection of 
the trial function. In practice, a fixed node approximation is 
used for fermions, such as electrons, which tests have shown 
give relatively small error when the nodes are obtained from 
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high-quality DFT orbitals for electronically simple materials 
such as c-BN. DMC for solids has been demonstrated to give 
significantly more accurate cohesive energies! 8], equations 
of state and Raman frequencies |9], and phase transitions! 10] 
than DFT. We have used both the CASINO QMC software 
suiteiH and QMCPACkIH. 

QMC simulations of solids are currently performed within 
the pseudopotential (PP) approximation, in which the core 
electrons are eliminated and their effect replaced by a nonlo- 
cal potential operator|8]. Since PPs are presently constructed 
with a lower-accuracy theory, such as Hartree-Fock (HF) or 
DFT, this replacement represents an uncontrolled approxima- 
tion. To eliminate this error, we develop a method for all- 
electron (AE) QMC simulations of soHds in QMCPACK using 
trial wave functions derived from full-potential linearized aug- 
mented plane wave (FP-LAPW) calculations using the EX- 
CITING code [13]. Space is divided into spherical muffin 
tin regions around the nuclei, and an interstitial region. Or- 
bitals are represented inside the muffin tins as a product of ra- 
dial functions and spherical harmonics, and outside as plane- 
waves. To ensure that a wavefunction satisfies the variational 
principle, it must be both continuous and smooth. We utilize a 
super-LPsP^ formalism that enforces continuity and smooth- 
ness at the muffin tin boundary. For efficiency, we represent 
the orbitals as 3D B-splines in the interstices and the product 
of radial splines and spherical harmonics in the muffin tins. 

Since AE QMC simulations are computationally expen- 
sive, we perform these simulations in 8-atom cubic supercells. 
Simulation cells this small would typically have significant 
finite- size errors. We eliminate these errors by combining data 
from AE simulations with that from PP simulations performed 
in both 8-atom and 64-atom supercells. Thus, we are able to 
simultaneously eliminate systematic errors from pseudpoten- 
tials and from finite- size effects. The corrected static-lattice 
energy is given at each volume as 

E = Ell + [El^ - El'] + A^/c + Ake'r^ (D 

where the term in brackets removes the psuedopotential 
bias. A^^^ and A|4^^^'^ are, respectively, potential! 14] and 
kinetic! 15] corrections for finite- size errors. We perform 
this procedure with three different PP sets commonly used 
in QMC: HF PPs from Trail and Needs (TN)!16, 17]; HF 
PPs from Burkatzki, Filippi, and Dolg (BED) [18]; and DET- 
GGA!^ PPs generated with OPIUM (WC)[20]. Perform- 
ing the same procedure with 128-atom PP simulations and 
16-atom AE simulations yields statistically indistinguishable 
results, demonstrating that finite- size errors are converged at 
this size. Additional methodological details and the finite- size 
data can be found in EPAPS document No. #. 

To compute the phonon free energy, we use density 
functional perturbation theory (DEPT) in the QUANTUM 
ESPRESSO package[21] with the Wu-Cohen functional and 
the OPIUM PPs. The phonon density of states, from which 
we derive thermodynamic data, is usually very well-described 
with DFT. 



(a) Parameters for 300K EOS 



Source 


PP/AE atoms 


Vo (A^) 


Bo (GPa) 


Q (nondim.) 


Trail-Needs PP + AE 
OPIUM GGA PP + AE 
BED HE PP + AE 


64/8 
64/8 
64/8 


11.792(18) 
11.769(17) 
11.781(20) 


381(6) 
385(6) 
382(7) 


3.87(6) 
3.86(6) 
3.87(7) 


BED HE PP + AE 


128/16 


11.812(8) 


378(3) 


3.87(3) 


Datchi et al.\24] 
Goncharov et al. 




11.8124 
11.817(32) 


395(2) 
387(4) 


3.62(5) 
3.06(15) 



(b) Parameters for thermal pressure 







^0(KA-^-) 


a^K-iA-^-) 


^n(KA-^") 


n = 





4.836656 xIQS 


2.598608x10-3 


-4.869563 x 10^ 


n = 


1 


-6.929704x101 


-1.200504x10-4 


1.763443 xl02 


n = 


2 


4.634278x10-1 


2.789128x10-6 


-2.186053 xlOO 


n = 


3 


-1.273468x10-3 


-2.008994x10-8 


9.303181 xlO-' 



(c) Parameters for Raman calibration 

co(cm-i) ci(cm-i) C2(K) h i?o(GPa) i?i(GPa) i?2(GPa) 
1055.9 -144.24 1497.8 3.0155 349.87 1849.4 112.33 



TABLE I: Parameters for the c-BN EOS and Raman calibration. 

We compute the free energy for our c-BN system at twelve 
unit-cell volumes, spanning volume compression ratios from 
0.84 to 2.0, corresponding to pressures of about -50 GPa to 
900 GPa. While it is not fully certain that the cubic phase of 
BN is stable to 900 GPa, no other structure has been observed, 
and theoretical studies have not identified any transitions be- 
low 1 TPa[22]. We use the Vinet formf^S*] for the isother- 
mal EOS, which we find represents our free-energy data very 
well, yielding the bulk modulus, Bq, its pressure derivative, 
Bq, and the equilibirum volume, (Table J[a)|). The statisti- 
cal error bars for each data point are directly determined from 
our QMC data. We compute statistical confidence ranges, tak- 
ing into account parameter cross-correlations with a simple 
Monte Carlo procedure. 

Figure [T] shows the EOS of c-BN at 300 K, with experi- 
mental data from Datchi et al [24] and Goncharov et al [7], as 
well as the present work from simulations with three different 
PPs. The residuals in (c) are derived from DMC simulation 
with PPs alone, while those in (b) combine all-electron and 
PP data. There is significant discrepancy between the theo- 
retical curves in (c), suggesting that PP simulation alone does 
not provide sufficient accuracy. Once the PP data is combined 
with AE data, as in (b), all the theoretical curves come into 
good agreement. Our theoretical EOS agrees reasonably with 
that in Refs. f27] andfT] within the experimentally measured 
pressure range, but the experimental extrapolation shows sig- 
nificant deviation at high pressure. 

We may write the thermal equation of state in the form 

PiV.T) = Psook{V) + Pth(^,T) - Pth(^,T = 300), (2) 

where P300K is the room-temperature contribution, fit to the 
Vinet form. The phonon contribution is written in an aug- 
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FIG. 1 : The c-BN EOS at 300 K, computed with DMC, compared to experiment. We plot the full EOS and the pressure residuals with respect 
to DMC with the BED PP. (a) gives the corrected EOS resulting from combining the EOS and PP data, while (b) gives the pressure residuals 
for the same data, (c) gives the uncorrected pressure residuals from the PP simulations only. Shaded areas represent the one-cr statistical 
confidence region from QMC data. 



merited Debye model, 

P,(VT) - J_E£^^ .3) 

e{V,T) = 0°{V)+^{V)exp{-a{V)T) (4) 
x{V) = xo+xiV+X2V'^ + X3V^, xe{e°,a,j3}{5) 

in which the Debye temperature, 9, is a function of both V 
and T (Table Ub)] ). The Debye free energy per two-atom cell, 
excluding the zero-point term, is given by 



Fd = QkeT 



In 



1 



-dx 



(6) 



c-BN can be used to calibrate pressure optically by mea- 
suring the frequency shift of the TO Raman mode, allowing 
bench-top experiments. We compute the pressure and temper- 
ature dependence of this frequency. Within the quasiharmonic 
approximation, phonon frequencies have explicit dependence 
on volume only. At constant pressure, however, these frequen- 
cies have an implicit temperature dependence resulting from 
thermal expansion. The dependence due to thermal expan- 
sion accounts for only about half the total T-dependence of 
the Raman mode, as observed in [25] and ||26|] . The remain- 
ing T-dependence can be attributed to significant anharmonic 
effects in c-BN, and is included in our calculations. 

Since the optical branch has small dispersion, we treat the 
anharmonicity as a one-dimensional on-site anharmonic oscil- 
lator, in a similar approach to the QMC computation of the TO 
Raman frequency of diamondJsl]. At each volume, we com- 
pute the effective Born-Oppenheimer potential well for the 
TO mode with DMC and the BFD PP at nine displacements 
along the mode eigenvector in the 64-atom supercell. We fit 
the data to a quartic polynomial, and numerically solve the ID 
Schrodinger equation in this analytic potential. This results in 
a set of single-phonon energy levels, {En}, with nonuniform 
separation. From {En}, we compute an intensity- weighted 
average Raman frequency, as a function of pressure and 



temperature. The matrix element for the transition from n 
to n — 1 is proportional to ^/n and is thermally weighted by 
the Boltzmann occupation of state n, so the intensity-averaged 
frequency, (u), can be given as 



he 



Z^n=l 



- , where In 



n e 



(7) 



The excess thermal softening, i.e. beyond that from thermal 
expansion alone, is accounted for by the thermal average of 
the anharmonic frequencies. Figure [2] shows the computed 
Raman frequencies compared with the experimental data re- 
ported in Refs. 125142711. 

Both Refs. ll25ll and ^2d\ give a ruby-like calibration for- 
mula, which can be expressed as 



{R/h) [{v/vf - 1] , 



(8) 



where R, v, and h have quadratic T-dependence. This de- 
pendence is sufficient below 2000 K, but cannot represent our 
data at high temperature. We use a form which captures the 
Boltzmann occupation of phonon excitations. 



v{P,T) = ^o(P) + ^i(P)exp 

nl/6 



MP) 



MP) 



bP 

Rn 



1 



n = 1,2,3 



(9) 



(10) 



with parameters in Table U^c)] and plotted in Figure [2(a)| Note 
that this formula cannot be analytically inverted, but a very 
simple iterative solution can be used for calibrating pressure 
from u and T. 

The main axis of Figure |2(b)| gives the room-temperature 
Raman frequency versus pressure. There is good agreement 
in the relatively low-pressure region in which the Raman fre- 
quency was measured. At very high pressure, the deviation 
with respect to the extrapolation in [25] increases with a max- 
imum discrepancy of 38 cm~^ or, conversely, a deviation in 
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FIG. 2: TO Raman frequency of c-BN as a function of temperature 
and pressure, (a) gives compares our theoretical Raman frequency 
with the fitted form in Eqs. [9]and[T0l (b) compares our caUbration 
at 300K with the experimental results from |25] and |26]. The black 
dots with error bars are the QMC frequencies, while the blue solid 
line gives the fit. The lower inset gives an expanded view at low 
pressure, in which the shaded region gives the statistical confidence 
region. The upper inset gives the variation with temperature. 



the pressure calibration of 50 GPa at 900 GPa. The devia- 
tions with respect to ll26ll are 70 cm~^ and 120 GPa. The 
experimental parameters capture the correct qualitative high- 
pressure behavior up to 900 GPa, despite the fact that data was 
available only to 20 and 64 GPa, respectively. This suggests 
the form for the fit was well-chosen. 

We have presented a fully ab initio pressure calibration 
based on quantum Monte Carlo simulations, and have intro- 
duced a method for all-electron simulations of solids to elimi- 
nate bias from pseudopotentials. This method should be appli- 
cable to at least first-row solids, allowing increased accuracy 
in the study of other materials and providing a new benchmark 
for other methods. The only remaining systematic error in the 
static contribution to the EOS is from the fixed-node approxi- 
mation used in DMC. For simple materials such as c-BN this 
error should be quite small, and tends to cancel between dif- 
ferent volumes. Thus we believe the EOS is robust enough to 
be used directly in experiment as a primary pressure calibrant. 



and can be used to cross-calibrate scales based on other ma- 
terials. Since the accuracy of our methods should not depend 
on compression, our calibration can be used up to 900 GPa. 
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